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Abstract 

In computing the spectra of quantum mechanical systems one en- 
counters the Fourier transforms of time correlation functions, as given 
by the quantum regression theorem for systems described by master 
equations. Quantum state diffusion (QSD) gives a useful method of 
solving these problems by unraveling the master equation into stochas- 
tic trajectories; but there is no generally accepted definition of a time 
correlation function for a single QSD trajectory. In this paper we show 
how QSD can be used to calculate these spectra directly; by formally 
solving the equations which arise, we arrive at a natural definition 
for a two-time correlation function in QSD, which depends explicitly 
on both the stochastic noise of the particular trajectory and the time 
of measurement, and which agrees in the mean with the ensemble 
average definition of correlation functions. 
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1 Introduction 



In quantum optics, a common experimental situation involves a system in 
a cavity which is monitored by measuring the spectrum of output photons. 
In order to compare theory to experiment it is necessary to calculate this 
output spectrum. 

What is commonly found is that the intensity of the output spectrum is 
given by the Fourier transform of two-time correlation functions of system 
variables In the case of a system isolated from its environment (apart 
from the measurement process itself), these correlation functions are simply 
the expectation values of products of Heisenberg operators at different times, 
e.g., {q{t2)q{ti)). 

As it is difficult to completely isolate a system from the environment, this 
approach does not always succeed. In the case of a quantum open system 
interacting continually with an external reservoir, one instead describes the 
system by a master equation, which - within the Markov approximation - 
can be written in Lindblad form: 

p = -t[H,p] +J2(Ln.pLl - l-{LlL^,p}). (1) 

The Lm are a set of environment operators which give the collective effects 
of the environment, and we have taken ^ = 1. 

The density operator p gives the probability for the expected outcomes 
of measurements on the system. In this case, one can still calculate output 
spectra, and find Fourier transforms of time correlation functions; these are 
no longer the expectation values of products of Heisenberg operators in the 
system Hilbert space alone, but are more complicated objects whose form is 
given by the quantum regression theorem [§, |]. 

Unfortunately, for complicated systems equation ([l|) can be very difficult 
to solve either analytically or numerically. In that case, it is often advan- 
tageous to consider an unraveling of the master equation into individual 
quantum trajectories, each one represented by a single state at every mo- 
ment in time. One of these unraveling techniques is quantum state diffusion 
(QSD) 0. In this, the normalized vector representing the pure state of 
the system evolves according to the QSD equation: 
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+ Y.iL^-{L^)mdCm. (2) 

This is an Ito stochastic differential equation, in which the dC,m represent 
independent complex Wiener processes. These satisfy 

M{dU)=M{dUdCn)=0, M{dC^d^n) = 5mndt, (3) 

where M represents an ensemble average of the noise. QSD reproduces the 
master equation in the mean: 

M(|^)(^|)=p. (4) 

(This is what is meant by an unraveling of the master equation.) Expectation 
values for operators obey a similar relationship: 

(0), = Tr{Op} = M((0)^). (5) 

The use of QSD as a practical algorithm to solve master equations has 
been widely investigated 0. This includes calculations of output spectra 
in quantum optics 0. While it seems reasonable that there should be a 
relationship between output spectra and time correlation functions in QSD 
analogous to that in Schrodinger or master equation dynamics, this has been 
difficult to show, as there is no generally accepted definition of a time cor- 
relation function for a single QSD trajectory. Since QSD is framed in terms 
of a nonlinear stochastic differential equation for an evolving state, it is not 
obvious how to generalize this to products of operators at different times. 

Gisin has attempted to provide a Heisenberg picture for stochastically 
evolving operators. Sondermann p, using the same stochastic equations, 
has suggested another definition for time correlation functions. In both cases 
it is difficult to relate this to the usual QSD formalism. In this paper we 
attempt to find a definition which arises from the QSD equation itself. 

In section 2 we derive the output spectrum of a quantum mechanical 
system, and show its relationship to time correlation functions, as given by 
the quantum regression theorem. 

In section 3 we derive a quantum output spectrum using QSD, and show 
how it leads to a natural definition of a two-time correlation function for a 
single QSD trajectory. Unlike the correlation functions which arise in the 
master equation derivation, these QSD functions depend explicitly on the 
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measurement time and the noise. On average, however, these dependencies 
vanish, and the result agrees with the usual definition of a time correlation 
function. 

In section 4 we look briefly at alternative unravelings, and discuss the 
use of these techniques for the practical computation of spectra. Finally, in 
section 5 we summarize our results and draw conclusions. 

2 Spectra and time correlation functions 

Consider a rather idealized model of an experimental set-up for measuring 
the output spectrum of a quantum mechanical system. Suppose we have 
a system with Hamiltonian Hi and Hilbert space Tii weakly coupled to an 
output mode with Hamiltonian H2 = tob^b and Hilbert space 7^2- We assume 
weak coupling to minimize the perturbation of the system evolution. For 
simplicity we assume a linear interaction: Hj = eq{b+b^) where the operator 
q represents a physical quantity of the system, like position for example, b 
and b^ are the annihilation and creation operators for the output mode, and 
e is small. 

The total Hamiltonian of the system plus output mode is 

H = H1 + H2 + HJ (6) 

and operates on the combined Hilbert space Hi^H2- Let the system degrees 
of freedom also be coupled to an environment, described by a set of environ- 
ment operators Lm acting on Hi. The system plus output mode obeys the 
master equation (|I|) with Hamiltonian (^. 

If the output mode is initially in the ground state, so that the initial 
density matrix is of the form p = poo ® |0)(0|, we can approximate at later 
times 

P ~ Poo ® |0)(0| + poi ® |0)(1| + pio ® |1)(0| + pii ® |1)(1| + 0{e'), (7) 

where the pij are time- dependent operators on Hi, we have explicitly sepa- 
rated the system and output mode degrees of freedom and neglected all but 
the flrst excited state of the output mode. We then have poi , pio ~ e and 
Pii ~ e^, and can rewrite equation (|l|) in the form 

Poo = >Cpoo + 0(6^), 
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poi = CpQi + iepooq + iup^i + O(e^), 
pio = CpiQ-ieqpQQ-iupiQ + 0{e^), 
pii = Cpn-ieqpoi + iepioq + 0{e'^). (8) 

£ is the time evolution superoperator restricted to the system degrees of 
freedom: 

Lpij = -i[Hi,pij] + [LmpijLl - ^{Ll^L^, pij}j . (9) 

To lowest order, poo evolves according to a normal master equation, with no 
reference to the output mode at all. The other components pij represent the 
weak signal transmitted via the interaction with the output mode. 
We can solve these equations: 

Poo(t) = ^oPoo(O); (10) 

poi(t) = ^e /* Sl,{poo{t')q)e^^^'-''^dt' = p\,{t) (11) 

Pu{t) = -te f Sl,{qpoi{t'))dt' + h.c. (12) 
Jo 



t nt' 



Jo 



4(g4',((5*"poo(0))g))e*"(*'-*")(it"(it' + h.c 



where 5**^ = exp(£(t2 — ^i)) is the time evolution superoperator from time 
ti to time t2 given by @. The excitation of the output mode is what is 
measured, so we are interested in the expectation value {b^b) = TT{b^bp} = 
Tr{pii}. We find this by taking the trace of equation Equation (|^) 

preserves the trace, so in taking the trace of (|1^) the factor exp(£(t — t')) 
has no effect, and may be dropped. Thus, 

{b^b){t) = 1^ 1^ Tr|g4',((5*"poo(0))g)}e^"(*'-*")rft"rft' + c.c. (13) 
^ Ij' Tr|gS'i:,(poo(t")g)}e^"^*'"*"^fii"rfi' + cc. 



f r'e^-(*'-*")(g(t')g(t"))QRTrft"t^t' + c.c. 
JO Jo 
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where the last equahty defines the two-time correlation function {q{t')q{t"))QKr- 
Note that (g(t')g(t"))QRT is determined by the evolution operator Sq of the 
reduced mixed state poo; in accordance with the quantum regression theorem 
[0]. We see that the expected output has the form of a Fourier transform of 
a two-time correlation function; this is just like the result in classical physics 
[^]. Examining output modes at different frequencies u gives the spectrum 
of the system. 

3 Spectra and time correlations in QSD 

We can calculate the results of section 2 with QSD by using relation (|^). 
We solve the QSD equation for a state in the combined Hilbert space 
H = Hi 7^2, starting with an initial condition |\E') = 1-00)10). By averaging 
{b^b)ys,(^t) over many trajectories we reproduce (|1^). This was done in P for 
the case of second harmonic generation. 

How can we interpret (6^6)ii,(j) for a single trajectory? Let's examine a 
little more closely the evolution (0). We can separate the components of \E' 

\m) = \Mt)m + \Hm) + o{e'), (u) 

where we again neglect all excited states of the output mode above the first. 
The initial condition is 

|0o(O)) = |^o), |0i(O)) = O. (15) 
The QSD equation becomes a pair of coupled equations 

d\4)o) = -iHi\(l)o)dt + J2(^{Lln)<t>oLm - ^Ll^L^ - ^{Ll^)^^{Lm)^o^\(po)dt 

+ J2(Lm- (L^)^o) \^o)dU + 0(6^), (16) 



i) = -iHi\(j)i)dt + Y,(^{Ll^)^^Lm - ^Ll^L^ - ^{Ll^)^o{Lm)^,^\(pi)dt 
+ J2(Lm- {Lm)4>^ \(t>i)dU - iuj\(pi)dt - ieq\^o)dt + ©(e^^), (17) 
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Note that ( [T6| ) is identical to the usual QSD equation for the system alone to 
first order in e. Just as in section 2, this approximation gains an extra order 
of e in accuracy for free. Thus, the system degrees of freedom have the usual 
quantum state diffusion behavior, essentially uninfluenced by the interaction 
with the output mode. Exactly as in (H), we have an almost unperturbed 
system, with a weak signal transmitted to the outside world. 

Equation ([T7|) is interesting, in that all of the expectation values in this 
equation are calculated with respect to |0o), so this is a driven linear equation 
with time-dependent coefficients. Given the solution |0o(^)), we can find 
\(f>i(t)), at least in principle. 

We formally integrate (|16D to get 

\Mt2)) = m,^o)ii\Mti)), (18) 

where T{$^,^o)ll time-evolution operator from time ti to time ^2- It 

explicitly depends on the noise ^ and the initial state \4>o), since the QSD 
equation is nonlinear. Given this time-evolution operator, the solution to 
(0)is 

\Mt)) = |V(e,^o)kl</'o(t'))e-'"^*-*')rft'. (19) 
The output spectrum is then 

(6t6)^(t) = f f\iJo\T^UT^ln,qn''\^o)e'^'^'^ (20) 
Jo Jo 

where the ^ and ipo have been suppressed for conciseness. This expression 
resembles a Fourier transform of some kind of correlation function, just as 
was the case in section 2; but this correlation function is defined for a single 
QSD trajectory. 

We can bring this rather closer to the treatment in section 2 by introduc- 
ing the projector 

=|^)(^|, 

= Foo|0)(0| + Poi|0)(l| + Pio|l)(0| + Pii|l)(l|, (21) 

where these partial projectors Pij are operators on Hi. These are related to 
our earlier treatment by 

Poo = |0o)(0o|, Poi = |0o)(0i|, Pio = |0i)(0o|, Pii = |0i)(0i|, (22) 
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and to section 2 by 

p., = M(P,,). (23) 

These Pij are not themselves projectors; however, P^^ = Pqo + 0{e^), and 
therefore can be considered a projector to good approximation. 

We define a time evolution superoperator for these partial projectors 

s{^,^lJoytlP = mi^oyti p T^{^,^oYti. (24) 

This is related to the superoperator for the master equation by 

M{S{^, i^oYtl) Poo = Sll Poo = exp{C{t2 - h)) Poo- (25) 
We shall need the following solutions: 

Poo{t) = S{C,ijo)oPoo{0), (26) 

Poi(t) = Z6 S{C, ,lJoYAPo,{t')qy-'^'~'^dt' = PUt) (27) 
Jo 

In order to compute the output signal {b^b)^(^t) we make use of the equality: 

Pii(t) = PioPoi. (28) 
In this way, we obtain the following expression for the output signal: 

(6t6)^(t) = rr{PioPoi} (29) 

= dt'TrS^Si^, iJoyAPoo{t')qy /J dt"S{^, ^o)^"(Poo(t")9)}e^"(*'-*"), 

= 1^ dt'TT!^T^,qPoo{t')Tf dt"r/„Poo(t")?rt''^}e*"^*'"*"^- 

Note that the two central operators T^Yt^,, can not be replaced with T^"\ 

Accordingly, in analogy with ([13D we define a measurement- dependent 
two-time correlation function for QSD: 

C(02,t2;Oi,ti;t) =Tr|^(e,V^o)L(02no(t2)) >9(?,Vo)UPoo(ti)Oi)} (30) 

Note that this assumes ti < t2', for ti > t2 an analogous expression can be 
formed. This function C in QSD has the nice feature that it is the trace of a 
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product of two operators, each involving one of the times ti and t2, similar to 
correlation functions of classical stochastic processes (which are products of 
the random variables at different times [Q). However, it also has the strange 
feature that the "final" time t appears in this definition. This is why we term 
C a measurement-dependent correlation function, rather than a true corre- 
lation function. Note that in perfect analogy with the quantum regression 
theorem, C{02, t2]0i,ti; t) is determined by the evolution operator S{^, ipoYtl 
of the reduced system's pure state |0o(^))- 

It remains to establish the relation between the function C{02, ^2; Oi, ti, t) 
and the QRT correlation function introduced in section 2. For this purpose 
we note that Pu = |0i)(0i| provides an equivalent expression for the output 
signal: 

{b^b)^it) = ^I{P^,} 

= Tr |5(e, i^oYt' {qS{^, i^,)UP,,{t")q)) }e-(*'-*")dtW + c.c. (31) 

Using this relation (^) we can derive an equivalent form of the function C: 

C{d2,t2;d,,h;t) = 

Tr|5(e,^o)L(o25(e,V'o)i?(i^oo(ti)Oi))}, (32) 

Note that in section 2 it was possible to remove the time evolution super- 
operator S^,, as it did not affect the trace. This is not true in the single 
trajectory explained above; so the function C depends explicitly on 

time, as mentioned above. Nevertheless, taking the mean over the noise of 
C{C)2,h', Oiiti'jt), one recovers the correlation function derived for mixed 
states from the quantum regression theorem: 

M(C(02,t2;Oi,ti;t)) = (02(t2)Oi(ti))QRT. (33) 

Note that the dependence on the final time t vanishes in the mean. 

While the definition (^) arises naturally in this derivation, and has the 
correct average behavior, its dependence on t remains a puzzling and rather 
annoying feature. It is possible to make a different definition, closely related 



to that of (|30|) , which avoids this problem. Since the operator S{^,'ijjo] 



t2 



in the definition (P^) vanishes in the mean, we can define a true two-time 



correlation function: 



(02(t2)Oi(ti))QSD = Tr|o25(e,^o)^?(Poo(ti)Oi)|. (34) 



9 



In this way the "final time" t disappears and the correlation function is 
formally identical to the QRT case (but for vectors instead of matrices). 
Moreover, the correlation function then appears as a scalar product of two 
vectors: 

(02(t2)Ol(tl))QSD = {T^^d,Mtl)\02Mt2))- (35) 

One can consider this new definition to be an average over "future" noise, 
i.e., noise after t2. 

(02(t2)Oi(ti))QSD = M{C{02,t2■,0^,h■,t))^^t'),t l<t2<t'<t- (36) 

Certainly, this once again reproduces the QRT correlation function in the 
mean, and gives exactly the same {b^b)^(^t) for a single trajectory: 

M((62(t2)Ol(tl))QSD) = (02(t2)Ol(tl))QRT. (37) 

While this definition removes the dependence on the final time t, this 
two-time correlation function still has some interesting features, related to 
the fact that the QSD time-evolution superoperator ( p^ does not preserve 
the trace. In particular, we note that if the second operator O2 is the identity 
J, then 

(/(t2)0i(ti))QSD 7^ (Oi(ti)), (38) 

contrary to the case of the QRT. If Oi = / no such difficulty arises. Since 
unravelings of the master equation do not generically preserve the trace, this 
feature will arise in any similar derivation. This might be one argument for 
using an unraveling specifically chosen to preserve the trace, such as that 
used by Gisin in his discussion of a Heisenberg picture for QSD [0, |[. In any 
case, the desired relation does hold in the mean: 

M((/(t2)0i(ti))QSD) = M((Oi(ti))). (39) 

Fortunately, correlation functions of the above form would not arise in any 
physically reasonable measurement scheme, since I does not describe an in- 
teraction. 

The definition of the two-time correlation function (0), and its relation 
to the quantum regression theorem (^) are the main results of this article. 
The central line of arguments can be summarized as follows: 

pn = M(Pn) = M(PioPoi) 7^ M(Pio)M(Poi) = PioPoi (40) 
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4 Alternative unravelings and practical cal- 
culations 



While the definition (|34D arises naturally from the QSD analysis of this prob- 
lem, it is not the only possible definition of a time correlation function for 
a single trajectory; indeed, there are an infinite number of such definitions, 
corresponding to different unravelings [|lO| . 

One such alternative unraveling has already been proposed by Gisin as a 
possible definition of a time-correlation function for individual trajectories, 
as well as providing a sort of Heisenberg picture corresponding to QSD [0]. 
Here we will suggest yet another such alternative unraveling. 

In particular, rather than solving the Lindblad master equation itself 
and then solving for the output spectrum, one might instead begin with the 
expression (|I3|) and attempt to unravel the time-correlation function directly. 

What form would such an unraveling take? In QSD, one unravels the 
density operator evolution into many trajectories, each consisting of a single 
state, and with a mean p = M{\iIj){iIj\). This works because p is hermitian: 
p = p^ . For a time correlation function, this is no longer sufficient. The form 
given by the quantum regression theorem is 

{q{t2)q{ti))QKT = Tr!^qSl-Mt,)q)y (41) 

where Sl^ is the time-evolution superoperator defined by the master equa- 
tion (0); however, the "initial state" p{ti)q is not hermitian. Therefore the 
evolution from ti to t2 cannot be unraveled in terms of pure states. 

One can, however, consider a pair of vectors 1-0) and such that M{\il)) 
does reproduce the correct evolution. Such a diad equation is quite analogous 
to the usual QSD equation. 

One pair of coupled equations that do the job are 

m = [-^H + {L^)^L - h^L - ^(Lt)^(L)^) \^)dt +{L- {L)^md^, 

{dct>\ = (01 (iH + Lt(L)^ - h^L - \{L^)^{L).^yt + (0|(Lt - {L'^)^)dC.m 

It is not hard to show that 

M(d|^)(0|) = /:(|V^)(0|), (43) 
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so this has the correct evolution in the mean. If IV') = |0) then these coupled 
equations reduce to the ordinary QSD equation (Q). Note that in general, 
however, the normalization of and |(^) is not preserved. 

The technique for calculating time correlation functions is as follows. If 
the initial density matrix is a pure state, p(0) = |'0o)('0o|) then one begins 
with both states equal \tl)) = |0) = \i/jo), and evolves them to time ti ac- 
cording to (1^), which is equivalent to the QSD equation @. At time ti, 
multiply {4>{ti)\ {(/){ti)\q and continue to evolve the diad according to p^). 
At time t2, multiply |'i/'(^2)) ?|'0(^2)) and take the trace. The mean over 
many such trajectories equals the time-correlation function (|^). 

This pair of equations (|^) shares many properties in common with the 
QSD equation (j^), but unlike QSD is not uniquely defined. Since only the 
composite diad |'0)((/>| is important, the norm and phase can be shifted arbi- 
trarily between these two states. Diosi has likened this to a gauge freedom 
[ p!0[| . Many such pairs of equations are therefore possible, as well as others 
with properties radically different from QSD; just as QSD is one of many 
unravelings of the master equation, albeit with unique symmetry properties. 
To solve for output spectra, of course, one must still Fourier transform the 
calculated time correlation functions. 

In fact, one can see that there are two distinct approaches to computing 
output spectra using QSD. One is to use the definition (Q) or alternative 
definitions such as (|4^) to calculate the time-correlation function, averaging 
over many runs, and taking the Fourier transform (|13|) . 

Alternatively, one can solve the QSD equation for the entire system plus 
output mode; the state is then in the larger Hilbert space Hi ® 7^2- This 
is essentially the approach taken by Schack et al. 0], who have also shown 
that these techniques can be used to calculate other quantities of interest, 
such as the spectrum of squeezing. 

Both approaches appear to have their advantages and disadvantages, and 
to be roughly equal in computational difficulty. It is likely that the best 
approach will vary from problem to problem. 

Note also that if we had unraveled the master equation (|I]) using an unrav- 
eling other than QSD, an exactly analogous argument would have followed. 
We would be able to resolve the equation for the system plus output mode 
into a pair of coupled equations, one corresponding to the unperturbed evo- 
lution of the system alone, the other to an output signal completely driven 
by the system. Thus, this type of argument could be used to define a notion 
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of time correlation functions for single trajectories in any unraveling, such as 



the Quantum Monte Carlo techniques or the orthojumps of Diosi 0, |TT|, |12 

These Quantum Monte Carlo (or Quantum Jump) techniques deserve 
further comment, as there have already been a number of papers published 
on their use in the calculation of time correlation functions |jl3|, |14|, p!5|] . 
These are relevant to our current discussion, since it can be shown that the 
equations for Quantum Jumps become identical to those for QSD in the case 
of heterodyne measurements |16[ 



The treatment of Gardiner, Parkins, and Zoller is particularly in- 
teresting in this context. Their definition of the time correlation function 
involves defining an auxiliary vector |/3,t) which is driven by the evolution 
of the quantum jump vector \(f),t), but does not in turn affect it. 

As stated in that paper [|1^ their equations are rather different. In par- 
ticular, they are considering only real noise, and have not gone to the weak 
coupling limit we have assumed in this paper, where the output has a neg- 
ligibly weak effect on the system over short times. Wiseman and Milburn 
[0 have generalized this treatment to consider the case of heterodyne mea- 
surements, introducing complex noise and showing that this limit is exactly 
equivalent to the Quantum State Diffusion equation. In the limit of hetero- 
dyne measurement with weak coupling to the external mode, the equation 
for 10, t) becomes the QSD equation (|^), and the two vectors t) and \(3, t) 
obey a pair of coupled equations identical to ([T6| ) and (|1^). The output spec- 
trum is given by the mean of t\(3, t), just as in (|20|) . From this, one could 
follow an argument exactly analogous to that of section 3 in this paper to 
arrive at a definition of a two-time correlation function identical to (^If). 



5 Conclusions 

In calculating quantum optical spectra, a common approach is to calculate 
the quantum time correlation function and derive the spectrum by taking its 
Fourier transform. This time correlation function has a form given by the 
quantum regression theorem, and requires a solution of the master equation. 

Quantum state diffusion provides in many cases an efficient method of 
solving the master equation. But hitherto, an appropriate definition of the 
time correlation function for a single QSD trajectory has been lacking, mak- 
ing QSD less useful for the calculation of spectra. 
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In this paper we have derived such a definition in a straightforward way, 
quite analogous to the derivation in the case of the full master equation. 
This correlation function in the mean has the form given by the quantum 
regression theorem, and the equations for it are very close to the original 
QSD equation, up to O(e^) in the interaction strength. 

This definition can be used as a practical numerical tool in computing 
quantum optical spectra. Other possible definitions and their potential for 
practical use have been briefly discussed. 

Finally, let us discuss the meaning of our result in the simple case of a 
damped harmonic oscillator at zero temperature. The stationary solution is 
the ground state. Hence, once the system has reached this state, nothing 
happens, and it is clear that the spectrum of the damped oscillator is not 
contained in the evolution of its state vector. Nevertheless, the quantum 
regression theorem tells us that all spectra, in particular the one correspond- 
ing to the position fluctuation, are contained in the evolution operator for 
the corresponding master equation. Similarly, the results presented here tell 
us that the spectra are also contained in the stochastic evolution operator 
^o)ti of the QSD description of the damped oscillator. The physics 
behind this is that whenever a spectrum is measured, the system's environ- 
ment is changed, hence its dynamics is perturbed. For example, to measure 
the spectrum of position fluctuations, something like weak position measure- 
ments have to be applied, and the ground state is no longer stationary [p!7| 



However, in contrast to standard quantum measurements, this perturbation 
can be made arbitrarily small (corresponding to small amplitudes of the mea- 
sured spectrum) over an arbitrarily long period of time. To flrst order, the 
system's evolution is unaffected, but its states acts like a source for the signal, 
as reflected by our equations (|16]) and (0). Hence the unaffected evolution 
operator (Liouville operator or QSD propagator or other stochastic propaga- 
tor 0, 0) contains the information about spectra that cannot actually be 
measured without affecting (i.e., weakly perturbing) the system's evolution. 
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Figure 1. The interaction between a quantum system with Hamiltonian 

Hi interacting with an external field mode H2 = ujh'^b via an interaction po- 
tential Hj = eq{b^ + b). The coupling of the system to the environment is 
modeled by a single environment operator L. Dotted lines represent Hamil- 
tonian terms and dashed lines represent coupling to the environment. 
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